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ABSTRACT 

An interesting new high-energy pulsar sub-population is emerging following 
early discoveries of gamma-ray millisecond pulsars (MSPs) by the Fermi Large 
Area Telescope (LAT). We present results from 3D emission modeling, including 
the Special Relativistic effects of aberration and time-of-flight delays and also 
rotational sweepback of B-field lines, in the geometric context of polar cap (PC), 
outer gap (OG), and two-pole caustic (TPC) pulsar models. In contrast to the 
general belief that these very old, rapidly-rotating neutron stars (NSs) should 
have largely pair-starved magnetospheres due to the absence of significant pair 
production, we find that most of the light curves are best fit by TPC and OG 
models, which indicates the presence of narrow accelerating gaps limited by ro- 
bust pair production - even in these pulsars with very low spin-down luminosities. 
The gamma-ray pulse shapes and relative phase lags with respect to the radio 
pulses point to high-altitude emission being dominant for all geometries. We also 
find exclusive differentiation of the current gamma-ray MSP population into two 
MSP sub-classes: light curve shapes and lags across wavebands impose either 
pair-starved PC (PSPC) or TPC / OG-type geometries. In the first case, the 
radio pulse has a small lag with respect to the single gamma-ray pulse, while 
the (first) gamma-ray peak usually trails the radio by a large phase offset in the 
latter case. Finally, we find that the flux correction factor as a function of mag- 
netic inclination and observer angles is typically of order unity for all models. 
Our calculation of light curves and flux correction factor for the case of MSPs 
is therefore complementary to the 'ATLAS paper" of Watters et al. for younger 
pulsars. 

Astrophysics Science Division, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA 

2 Unit for Space Physics, North- West University, Potchefstroom Campus, Private Bag X6001, Potchef- 
stroom 2520, South Africa 

3 NASA Postdoctoral Program Fellow 

4 Max-Planck-Institut fur Radioastronomie, Auf dem Hugel 69, 53121 Bonn, Germany 



-2- 



Subject headings: acceleration of particles — gamma rays: theory — pulsars: 
general — radiation mechanisms: non-thermal — stars: neutron 



Introduction 



The field of gamma-ray pulsars has already benefited profoundly from discoveries made 
during the first year of operation of the Fermi / Large Area Telescope (LAT). These in- 
clude detections o f the radio-quiet gamma-ray pulsar inside the supern ova remnant CTA 1 



( lAbdo et al.l 120081 ) , the second gamma-ray millisecond pu 



following the EGRET 4.9a- detection of PSR 



confidence EGRET pulsars ([Thompson et al. 
radio-quiet pul sars using b 



been unveiled (Abdo et al. 



ind searches (Abdo et al. 



0218 



1999 



sar (MSP) (jAbdo et al.l l2009ah 



4232 (IKuiper et al.ll2004h . the 6 high- 



Thompsonl I2004T ). and discovery of 16 



2009bf ). In addition, 8 MSPs have now 



2009dl. see Table [D), confirming expectations prior to Fermi's 



launch in June 2008 (IHarding. Usov. fc Muslimovil2005l ; IVenter fc De Jagerll2005br) . A Fermi 
six-month pulsar catalog is expected to be released shortly (jAbdo et al.ll2009el ). AGILE has 
also reported the discovery of 4 new gamm a-ray pulsars, and marginal detection of 4 more 



(IHalpern et al 



pulsars ( 



2008 



Pellizzoni et al. 



Pellizzoni et al.ll2009bf ). in addition to the detection of 4 of the EGRET 



2009al ). Except for the detection of the Crab at energies above 

jeen detected by 



25 GeV (lAliu et al.ll2008h . no other pulsed emission from pulsars has as yet 



ground-based Cherenkov telescopes (ISchmidt et al. 



2005 



Albert et al 



2007l:lFussling et al.ll2008l : iKildealboQsl : iKonopelkjliooi Albert et al. 
De los Revesll2009h . 



2007 



Aharonian et al. 



20081 : ICelik et al. 



200. 



MSPs are characterized by relatively short periods P < 30 ms and low surface magnetic 
fields B ~ 10 8 — 10 9 G, and appear in the lower left corner of the PP-diagram (with P the 
time-derivative of P; see Figure HJ where the newly-discovered Fermi MSPs are indicated by 
squares). MSPs are thought to have been spun- up to millisecond periods by tra nsfer of mass 



and angular momentum from a binary companion during an accretion phase (lAlpar et al. 
19821 ) . This follows an evolutionary phase of cessation of radio emission from their ma- 



ture pulsar progenitors, after these have spun down to lon g periods and crossed the "d eath 
line" for radio emission. These "radio-silent" progenitors ( jGlendenning fc Weber! 120001 ) are 
thought to reside in the "death valley" o f the PP-diagram, which l ies below the inverse 
Compton scattering (ICS) pair death line (IHarding fc Muslimovll2002l ). 



The standard "recycling scenario" (IBhattacharya &: Van den Heuvellll99ll ) hypothesiz- 
ing that MSP birth is connected to low-mass X-ray binaries (LMXRBs) might have been 
confirmed recently by the detection of ra dio pulsations from a nearby MSP in an LMXRB 
system, with an optical companion star (jArchibald et al.l 120091 ). Optical observations indi- 
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cate the presence of an accretion disk within the past decade, but none today, raising the 
possibility that the radio MSP has "turned on" after termination of recent accretion activity, 
thus providing a link between LMXRBs and the birth of radio MSPs. 



two emission regions. Polar cap (PC) models ( 


Harding et al. 


1978: 


Daughertv & Harding 


1982; 


Sturner et al. 


1995; 


Daughertv & Harding Il996 


) assume extraction of primaries from 



the stellar surface and magnetic pair production of ensuing HE curvature radiation (CR) 
or ICS gamma rays, leading to low- altitude pair formation f ronts (PFFs) which screen 



the accelerating electric field (IHarding &: Muslimovlll998l . l200ll 120021 ) . These space-charge- 



limited-flow (SCLF) models have since been extended to allow for the variation of the CR 
PFF altitude across the PC and therefore accel eration of primaries along the last open 



magnetic field lines in a slot gap (SG) scenario (lArons fc Scharlemannl Il979l ; lArond 11983 
Muslimov fc Harding 2003 . 2004a ; Harding et al. 20081 ). The SG results from the absence of 
pair creation along these field lines, forming a narrow acceleration gap that extends from the 
neutron star (NS) surface to near the light cylinder. Th e SG model is thus a possible physi- 
cal realization of the two-pole caustic (TPC) geometry (iDyks &; Rudakl 120031 ) . developed to 



explai n pulsar HE light curves. On the other hand, outer gap (OG) models (ICheng et a 
1986al lbl: Ichiang fc Romanilll992l . 1 19941 : lRomani|[l~9~9~6l ; ICheng et allboool : Izhang et al.ll2004h 
assume that HE radiation is produced by photon-photon pair production-induced cascades 
along the last open field lines above the null -charge surfaces (O • B = 0, with Q = 2tt/P), 
where the Goldreich- Julian charge density (jGoldreich &; Julian! 119691 ) changes sign. The 



pairs screen the ac celerating E-field, and limit both the parallel and transverse gap size 
(ITakata et al.ll2004l ). Classical OG models may be categorized as "one-pole caustic models", 
as the assumed geometry p revents observ ation of radiation from gap s (caustic s ) associated 
with both magnetic poles (iHardingi 120051 ) . Mo re recently, how ever, iHirotanil (120061 . 120071 ) 
found and applied a 2D, and subsequently a 3D (jHirotanill2008bl ) OG solution which extends 
toward the NS surface, where a small acceleration field extracts ions from the stellar sur- 



facein an SCLF-regime (see also ITakata et al.l (120 04. 20Q 6|), an d in particular ITakata et al. 



J2008I ) for application to Vela). Lastly, ITakata & Chanel J2009h modeled Geminga using an 
OG residing between a "critical" B-field line (perpendicular to the rotational axis at the 
light cylinder) and the last open field line. 

Current models using dipole field structure to model MSPs predict largely unscreened 
magnetospheres due to the relatively low B-fields inhibiting copious magnetic pair produc- 
tion. Such pulsars may be described by a variation of the PC model (appl icable for younger 
pulsar s ), which we will refer to as a "pair-starved polar cap" (PSPC) mo del (IMuslimov fc Harding 
2004bl ; IHarding. Usov. fc Muslimovi I2OO5I : IMuslimov & Harding! [20091 ) . In a PSPC model, 
the pair multiplicity is not high enough to screen the accelerating electric field, and charges 
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are continually accelerated up to high altitudes over the full open-field-line region. The 
formation of a PSPC "gap" is furthermore naturally understood in the context of an SG 
accelerator progressively increasing in size with pulsar age, which, in the limit of no electric 
field screening, relaxes to a PSPC structure. 



Sever al authors have modeled MSP gamma-ray fluxes, spectra and lig 



the PS PC dFrackowiak fc RudajJ2005 aB lHarding. Usov. fc Muslimov 



2005a: Venter 



2008 



it curves in both 



Zaiczvkil2008h and OG flZhang fc Chene 



2003 



Colle c tive emission from a population of MSPs in globular clusters ([Harding. Usov. fc Muslimov 



2005 



Zha ng et al. 



Venter &; De Jager 
2007h 



cases. 



2005 



Zhang et all 120071 : Bednarek k Sitarekl 120071 : IVenter fc De Jager 



2009J) and in the Galactic Center (IWangj 120061 ) have also been considered. 



2008 



Venter et al. 



Watters et al. 



(120091 ) recently calculated beaming patterns and light curves from a population of canonical 



pulsars with spin-down luminosities E TOt > 10 ergs - using geometric PC, TPC, and OG 
models. They obtained predictions of peak multiplicity, peak separation, and flux correction 
factor fa as functions of magnetic inclination and observer angles a and (, and gap width w. 
The latter factor fa is used for converting observed phase-averaged energy flux G b s to the 
total radiated (gamma-ray) luminosity L 7 , which is important for calculating the efficiency 
of converting E rot into L T A good example is the inference of the conversion efficiencies of 
globular-cluster MSPs which ma y be collectively res ponsible for the HE radiation observed 
from 47 Tucanae by Fermi-LAT ( lAbdo et al.ll2009d ). 



In this paper, we present results from 3D emission modeling, including Special Relativis- 
tic (SR) effects of aberration and time-of-flight delays, and rotational sweepback of B-field 
lines, in the geometric context of OG, T PC, and PSPC puls ar models. We study the newly- 
discovered gamma-ray MSP population (lAbdo et al.l l2009dj ) . and obtain fits for gamma-ray 
and radio light curves. Our calculation of light curves and flux c orrection facto r s foia , (, P) 
for the case of MSPs is therefore complementary to the work of IWatters et al.l (120091 ) which 
focuses on younger pulsars, although our TPC and OG models include non-zero emission 
width. Section [2] deals with details of the various models we have applied. We discuss light 
curves from both observational and theoretical perspectives in Section [3X and present our 
results and conclusions in Sections 0] and [51 



2. Model Description 



2.1. B-field and SR Effects 



Deutschl (119551 ) found the solution of the B- and E-fields exterior to a perfectly-conducting 



sphere which rotates in vacuum as an inclined rotator. We assume that this retarded 
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vacuum dipolar B- field is representati ve of the magnetospheric structure, and we use the 
implementation by iDvks et al. (l2004aT); iDvks Sz Harding (2004) , following earlie r work by 



Romani &: Yadigaroglu 



(1995) 



Higgins & Henriksenl(ll997n : lArendt fc Eilekl (ll998f ): ICheng et al 



( 120001 ). For this B- field, the PC shape is distorted asymmetrically by rotational sweepback 



of field lines. Each field line's foo tpoint is labeled by the open volume coordinates (r ovc , l 



as defined by IDvks et al.l (j2004al ) , with r ovc labeling self-similar contours or "rings" (r ovc is 
normalized to the PC radius Rpc), and / nvn giving the arclength along a ring (analogous to 
azimuthal angle; also refer to lHarding et al.l (120081 ) for more details). 



We calculate the rim of the PC by tracing field lines which close at the li ght cylinder back 
to th e stellar surface, and then divide this PC into rings (see e.g. Figure 2 of lDyks &: Harding 
20041 ) and azimuthal bins, with each surface patch dS associated with a particular B-fleld 
line. We follow primary electrons moving along each field line, and collect radiation (cor - 
rected for SR-effect s ) in a phaseplot map ( Section 12. 3j) . Following IChiang fc Romanil (119921 ); 
Cheng et al.l (120001 ); IDvks fc Rudakl (120031 ) . we assume constant emissivity along the B-lines 
in the gap regions of the geometric PC, OG, and TPC models (but not for the PSPC model), 
so that we do not need to include any particular E-field (or calculate dS explicitly) for these. 



In the case of the PSPC model, we use the approximation £ 



(with £ = 9/9 pc the 



normalized polar angle, and 8 pc 
to high altitudes (Section 12.21) . 



{ytR/cfl 2 the PC angle), and include the full E-field up 



In addition to the rotational sweepback (retardation) of the B-lines, we include the 
effects of aberration and time-of-flight delays. We calculate the position and direction of 
photon propagation (assumed to be initially tangent to the local B-line) in the co-rotating 
frame, and then aberrate this direction using a Lorentz transformation, transforming from 
the instantaneously co-moving frame to the IOF. Lastly, we correct the phase at which the 



photon reaches the observer for time delays due to t he finite speed of 



about calculation of these SR effects ma y be found in lDvks et al.l (l2004bl); 



( 20041 ). following previous work by e.g., iMorinil ( 19831 ); iRomani &; Yadigarog kJ Jl995h . We 



it. More details 



Dvks fc Harding 



furthermore explicitly use the curvature radius of the B-field lines as calculated in the inertial 
observer frame (IOF), and not in the co- rotating frame, when performing particle transport 
calculations (Section I2.2p . Such a model has also recently been applied to the Crab by 
Harding et all (b()08h . 



We have calculated TPC and OG models assuming gaps that are confined between 
two B-field lines with footpoints at r ovc l and r ovc 2 . We therefore activated only a small 
number of rings near the rim (r ovc ~ 1) with r ovc G [tw^, r OVCi 2], and binning radiation 
from these, assuming constant emissivity over the emitting volume. For TPC models, we 
used r ovc e [0.80,1.00], [0.60,1.00], [0.90,1.00], [0.95,1.00], and [1.00,1.00] (see Table [2]) 
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corresponding to gap widths of w = r OVCi2 — r OVCj i = 0.20, 0.40, 0.10, 0.05, and 0.00 (the 
last of these is what we have referred to as the TPC model). Similarly, we investigated 
OG models with r ovc e [0.90,0.90], [1.00,1.00], [0.95,1.00] (widths of w = 0.00, . 00, ari d 
0.05). These widths are smaller than e.g. the value of ~ 0.14 used by IHirotanil (j2008al ). 
We did not find good light curve fits for TPC models with large w. In the case of OG 
models, one should consider non- uniform emission when choosing large w, which is beyond 
the scope of this paper. The assumption of constant emissivity in the emitting volume is 
a simplification, as OG models are expected to produce the bulk of the gamma-radiation 



along the inner edge 



ovc, inner 



< r ovc < r OVCjPFF ) of the gap (r OVCjPFF < r ovc < 1), with 



r ovc,PFF indicating the position of the PFF, and r nv r ,j me . r some smaller radius depending on 
the radiation surface thickness (IWatters et al.ll2009l ). We lastly modeled the PC and PSPC 
cases with r ovc 6 [0.00, 1.00] (i.e., the full open-field-line volume, for both constant emissivity 
and full radiation codes). We used 180 colatitude (£) and phase (4>) bins and individual ring 
separations of 5r ovc = 0.005, while collecting all photons with energies above 100 MeV (in 
the case of the PSPC model) when producing phaseplots and subsequent light curves. 

It is important to note a critical difference bet ween the radiation distribution in our TPC 
and OG models and that of IWatters et al.l (120091 ). We assume that emission is distributed 
uniformly throughout the gaps between r ovc l and r OVC)2 , so the radiation originates from 
a volume with non-zero width across field lines and the radiation and gap widths are th e 
same. Fo r the TPC model, t his geometry is similar to that adopted by lDyks et al.l (j2004al ). 
although iDyks et al.l (j2004al ) assummed a Gaussian distribution of emission centered at the 
gap midpoint while we simply assume a constant emissivity a cross the gap, both of which 
crudely approximate the radiation pattern expected in the SG. IWatters et al.l (120091 ) assume 
that the emission occurs only along the inner edge of both the TPC and OG gaps (r ovc ,i in 
our notation), and so their radiation width is confined to a single field line and not equal to 
their gap width (w in their notation). In the case of the OG, the physically realistic emission 
pattern would h aye a non-zero w idth lying somewhere between infinitely thin and uniform 
assumptions (see IHirotanil l2008bl ). 



2.2. Particle Transport and PSPC E-fleld 

We only consider CR losses suffered by electron primaries moving along the B-field lines 
when modelin g the HE emission. In this case, the (sin gle electron) transport equation is 
given by (e.g., ISturnerlll995l ; iDaugherty & Harding! 1 19961 ) 



E P 



-^e,gain 'TcR^eC 



e(3 r cE\\ 



2fc 



(1) 
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where c is the speed of light in vacuum, f3 r = v r /c ~ 1 the particle velocity, e is the electron 
charge, 7 is the e lectron Lorentz f actor, 7 CR m e c 2 the frequency-integrated (total) CR loss 
rate per particle (IBulik et al.ll2000i ). p c the curvature radius (as calculated in the IOF; see 
Section f2.1|) . and E\\ the accelerating E- field parallel to the B-field. T he acceleration and 
loss terms balance at a particular 7rr in the radiation reaction regime (ILuo et al.ll2000l ): 



7rr 



3^ 
2e® 



(2) 



Previous studies (e.g., Venter fc De Jager 2005a ; Frackowiak fc Rudakll2005a b; Harding, Usov. fc Mus 



2005 



Venter fc De Jagerj|2008 ; jZaiczv kl l2008h have used the solutions of lMuslimov fc Hardin 



(Il997h : lHarding fc Muslimovi fll998h for the PSPC E-field: 



(i) 



E, 



(2) 



--^ (^ R ) {12/c'si cos a + 6s 2 8$ R H(l)6'(l) sinacos0 pc } 



£0 
R 



W R ) 2 



3 k' 
2^f 



cos a + ^^(^^(^^(^esinacos^pej (l - £ 2 ) , 



(3) 

(4) 
(5) 



with 



$ 



e = 



e GR ( v ) 



Si 



S2 



1/2 



'pc 



B^VtR? 

5 

c 

2GM 
ttR T] 



^ GR (r/)(l 

A"; 



0^1/2 



77^ GR (?7)(1 - e/77)V2 

1 _ e -7(l)(»J-l) 



(6) 
(7) 



(9) 
(10) 

(11) 

(12) 
(13) 



and ki and ki are the positive roots of the Bessel functions Jo and J\ (with ki + \ > ki and 
k i+ i > ki)] 9q R = 6* GR (1); 7(1) may be 7»(1) or 7i(l) in the expression for T\. The functions 
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H(rj), f(rj), and S'(rj) are all of order unity, and are defined in iMuslimov fc Tsyganl (119921 ). 



The first solution is valid for rj - 1 < 1, and £, ( , 2) for 6^ R < 77 - 1 < c/(QR); R is 
the stellar radius, rj = r/R, a the angle between the rotation and magnetic axes, <ft pc the 
magnetic azimuthal angle, k' = 2GI/(c 2 R 3 ) the General Relativistic (GR) inertial frame- 
dragging factor (distinct from the k(x) function to be defined later), and / the moment of 
inertia. 



Muslimov fc Harding! (j2004bl ) found the solution of En for altitudes close to the light 



cylinder in the small-angle approximation (small a, £, and high altitude): 

3 



E 



(3) 



3 (ttRY B 



16 V c ; /(l) 



K | I 

1/2 



Vc 



cos a 



D(— J A(l + 2£ 2 ) 
x £sinacos0 pc ] (l — £ 2 ) , 



(14) 



and A is defined after Eq. (35) of IMuslimov fc Harding! (j2004bl ). They proposed that one 
should employ the following formula to match the last two solutions: 



En 



E^exp [—(t] - 1)/ (t] c — 1)]+E 



(3) 



(15) 



with r] r a radial parameter to be determined using a matching procedure. IMuslimov fc Harding 



(12004bf ) estimated that r] c ~ 3 - 4 for MSPs when £ 



GR 



0.5. 



(3) 

It is important to include the high-altitude solution Er, , as Fermi resu lts seem to 



indica te that the HE radiation is originating in the outer magnetosphere (e.g., lAbdo et al. 
2009dl ). Beaming properties and spectral characteristics of the emission may therefore be 



quite different in comparison to calculations which only employ £n and E^' . In addition, 
while we use E-field expressions derived in the small-angle approximation, it is preferable to 
use the full solution of the Poisson equation, particularly in the case of MSPs which have 
relatively small magnetospheres and therefore much larger PC angles compared to canonical 
pulsars. 

In this paper, we calculate rj c (P, P, a, £, pc ) explicitly for each B-field line according 
to the following criteria (we use P = 1(T 20 , M = 1.4M , R = 10 6 cm, and / = OAMR 2 
throughout). We require that the resulting E-field should: 

(1) Be negative for all 1 < rj < c/(flR); 

(2) (3) 

(2) Match the part of the E\ -solution which exceeds En in absolute magnitude (i.e., where 



n(2) 



—Eu > —E^') as closely as possible; 
(3) Tend toward E^ for large rj. 



n(3) 
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The first criterion is required to mitigate the problem of particle oscillations which occurs 
when the E-field reverses sign beyond some al titude. Instead of this happening ~ 40% of the 
time (jVenter fc De Jagerll2005al ; I Venter! 120081 ) . we now only have to ignore solutions where 
— E\\ < for i] > 1.1 for ~ 8% of the time. The two lower-altitude solutions ErP and 
have been matched at rj = rjb, using (jVenterll2008l ) 



r]b 



1 + 0.0123P" 0333 . 



(16) 



Example fits for E\\ are shown in Figure [2] for different parameters, as noted in the caption. 
The top two panels show fits for two different T] c , while the bottom panel is an example where 
no solution for r\ c is found (according to the first criterion above). 

For illustration, Figure [3] shows contour plots of r] c ~ 1 — 6 for different a, £, and 
4>p C , and for P = 5 ms; £ is the 'radial' and pc the azimuthal coordinate for these polar 
plots. From these plots, one may infer that the "oscillatory solutions" are encountered when 

(2) 

4> pc ~ 180°, and for large a (which is where the second term of E^ becomes negative and 
dominates the first positive term inside the square brackets of Eq. [1]). The resolutions 
become progressively smaller for these cases, until no solution is found which satisfies the 
above criteria; we ignore emission from those particular field lines. 

We tested our full solution of E\\, which incorporates E^ through E^\ for conservation 
of energy when solving the transport equation (Eq. [T]) for relativistic electron primaries. 
Fi gure HI indicates the log^Q of acceleration rate 7gain — E e ^ a j in /Tn e cP' ) loss rate 7i oss = 7 CR , 
curvature radius p c , and the Lorentz factor 7 as functions of distance. Although we did not 
find perfect radiation reac tion wh e re the acceleration and loss terms are equal in magnitude 
(similar to the findings of IVenterl (120081 )). integration of these terms along different B- field 
lines yielded energy balance (i.e., conversion of electric potential energy into gamma-radiation 
and particle kinetic energy) for each integration step of the particle trajectory. An example 
of this is shown in Figure [51 where the graph of the cumulative energy gain (f^ =1 ^7 ga in) 
coincides with that of the sum of the cumulative energy losses and the acquired particle 
energy (/^ =1 <i7ioss + j(v) ~ 7o) f° r & U V (to within ~0.3%), with 7 = 7(77 = 1) the initial 
Lorentz factor at the stellar surface. We used 70 = 100, but the calculation is quite insensitve 
to this assumption, as 7 quickly reaches values of ~ 10 7 (FigureHJ). The quantities in Figure[5] 
are plotted in units of m e c 2 . 
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2.3. Generation of Phaseplots 



In the case of the PSPC model, we normalize the particle outflow along each B-line 
according to 

p e (f] = l,C,0pc) 



dN{£, 



'PC) 



-dSp c, 



(17) 



with dN the number of particles leaving a surface patch dS per unit time with initia l speed 
Pqc and p e is the charge density given by Eq. (12) of lHarding &: Muslimovl (119981 ). The 
latter is equal to the GR equivalent of the Goldreich- Julian charge density at the NS surface. 
The express i on in Eq. (|17|) is similar to the classical Goldreich- Julian expressions used by 
Storv et all (ho07h : 



N, 



G.J 



™GJ 



1.3 x 10 30 B 12 P- 
2ir (1 — cos 9 



particles s, 



(18) 
(19) 



with N Gi the total number of particles injected from the PC per unit time, B 12 = B /10 G, 
77-gj the injected particle flux, and dn^i = ncjdS analogous to the GR quantity dN defined 
in Eq. (|T7|) . While the classical injection rate dncj is constant across the PC, the GR 
expression we use has both a £- and pc -dependence. (Even though dN varies across the 
PC, we assume that it stays constant along B-lines, i.e. that it has no 77- dependence.) 

We have distributed primary electrons uniformly across the PC using a constant step 
length dl ovc along all rings between consecutive electron positions, so that there are generally 
less electrons per ring for the inner rings than for the outer ones. Because of this uniform 
distribution, we could approximate the area of the surface elements using 



dS 



iV e ,tot ' 



(20) 



with iV e tot the total number of electrons positioned on the PC surface (depending on grid 
size of the mesh into which the PC area was divided). These electron positions coincide with 
B-line footpoints on the stellar surface. We next followed the motion of electron primaries 
along these lines (Section l2.2p . collecting HE radiation and binning as described below. 



The instantaneous CR power spectrum is given by (lJacksorull975l ; lHardinpl981l ; iDaugherty fc Harding 
19821 : Istorv et~al1l2007h 



\dEj CR 



— ) = V^GJfin 



a 



with 



3A c7 3 



2-rrp, 



2p c 2m e c 2 p c 



(21) 



(22) 
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the c ritical energ y, A c = h/(m e c) the Compton wavelength, a^ ne the fine-structure constant, 
and (jEr perl 1 19661 ) 



= U 5/3(X } l 1-253 *W x » 1, (23) 

with the modified Bessel function of order 5/3. We calculate the number of CR photons 
radiated per unit time by the primaries in a spatial step cisioF (as measured along the B-field 
line in the IOF), in an energy bin of width dE = E2 — E\, using 

dn 7m = %^^iV, (24) 

with 

E - = w (Bl + s) ' (25) 

and 

f5 2 k(x) dx 

W = J * ; , , (26) 

with £?o <C 1. The expression in Eq. (12^1) gives the number of photons radiated per primary 
per unit time with an energy ~ i^bin (i-e., the ratio of power radiated per primary in a 
particular energy bin to average bin energy, in m e c 2 units) multiplied by a time step dsioF / c, 
multiplied by the number of primaries passing per unit time dN. (The 'weighting factor' 
W therefore scales the total power to the power radiated in the particular energy bin.) We 
ignore field lines with dN < 0. 

For all the other geometric models, we assume constant emissivity per unit length, i.e. 
dh 7> cR oc dsioF- 

We lastly accumulate <in 7i cR in (C> 0)-bins (after applying the SR effects described in 
Section |2~TTI) . and divide by the solid angle subtended by each phaseplot bin, dfl = (cos( — 
cos(( + d())dcf) ~ sm(d(d(f>, to make up the final phaseplot. 



2.4. Radio Beam Model 

We model the radio emission beam using an empirical cone model that has been de- 
veloped over the years through detailed study of pulse morphology and polarization char- 
acteristics of the average-pulse profile. The average-pulse profiles are quite stable over long 
timescales and typically show a variety of shapes, ranging from a single peak to as many as 
five separate peaks. The emission is also highly polarized, and displays changes in polariza- 
tion position angle across the profile that often matches the position angle swing expected 
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for a sweep across the open field lines near the magnetic poles in the Rotating Vector Model 
( IRadhakrishnan fc Cookdll969l ). 



Rankin's (lRankinlll993l ) study of pulse morphology concluded that pulsar radio emission 
can be characterized as having a core beam centered on the magnetic axis and one or more 
hollow cone beams also centered on the magnetic axis surrounding the core. Although 
Rankin's model assumes that emission fills the core and cone beams, other studies (e.g., 



Lyne fc Manchester! Il988l ) conclude that emission is patchy and only partially fills the core 
and cone beam patterns. 



The particular desc ription we adopt is fromlGonthier et al.l (120041 ) and is based on work 
of Arzoumanian et al. (lArzoumanian et al.l 120021 ) . who fit average-pulse profiles of a small 
collection of pulsars at 400 MHz to a core and single cone beam model based on the work 
of Rankin. The flux from the con al component seen at angle 6 to the magnetic field axis 
(modified by iGonthier et al.l (120041 ) to include frequency dependence v) is 



s(e,v) 



The annulus position and width of the cone beam are 

6 = (1.0 -2.63 S w ) Pcone , 



where S w = 0.18 ([Harding et al.l 120071 ). and 



n — 1 O/|° r 0.5 p-0.5 



with 



r KG ~ 40 



P 



0.07 



lO-^ss" 1 



p0.3 -0.26 
r GHz, 



(27) 

(28) 
(29) 

(30) 
(31) 



the radio emission altitude in units of stellar radius (IKijak fc Gilll2003l ). and z^ghz = f/1 GHz. 
(We do not assume a longitudinal extension of the radio emission region, but only use a single 
emission altitude.) According to Eq. (I3TI) . the altitude of the conal radio emission is a weak 
function of P, but the emission occurs increasingly close to the light cylinder (at -Rlc = c/Q) 
as P decreases (for more or less constant P) . For Crab-like periods, the conal emission occurs 
at altitudes of 10% — 20% of the light cylinder radius (and similar for typical MSP parameters 
of P ~ a few milliseconds and P ~ 10~ 20 ). For the current study, we are only interested in 
pulse shapes and phase shifts between the radio and gamma-ray pulses. We therefore use 
relative units for the cone beam luminosity. 
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2.5. Flux Correction Factor 



It is very important to be able to scale from the observed (phase-averaged) energy flux 
G bs to the all-sky luminosity, as this is used to define the gamma-ray radiati on efficiency 



7 7 7 , a cr ucial quantity in characterizing the energetics of pulsar emission (see e.g. lAbdo et al. 



( I2009el ). where 7/ 7 oc /q). Such a flux correction factor is necessarily model-dependent, 
as any observer only sees a small part of the total radiation: that coming from a slice through 
the emission beam, determined by the line-of-sight (. 



Venter! (120081 ) defined the total gamma-ray luminosity using 



Ad 2 G b s , 



(32) 



rbeam 



with A = eAtt /f3 ohs , e = /? obs G beam /G' obs , (3 ohs t he duty cycle Att" the average 
beaming angle, and G beam the all-sky total energy flux. I Watters et al.l (120091 ) used a similar 
definition 

(33) 



/n(a,Ci 



v 7 — 4:71 f q(1 2 G ODS , 

II F 1 (a,(,(j))sm(d(d(f) 



(34) 



2/F 7 (a,CM 

with F 7 the photon flux per solid angle ('intensity'), and (e the Earth line-of-sight, so that 



A ~ 4vr/ r , 



(35) 



assuming similar distributions of gamma-ray photon and energy fluxes in ((, 0)-space (i.e. 
F 7 (C, 0)/F 7)tot « G 7 (C, <fi)/G 7)tot ). In Section HI we calculate fa for different pulsar models, 
using Eq. 



3. Light Curve Data 



We compare the light curves generated with the different theoretical models (by making 
constant-^ cuts through the respective phaseplots of gamma-ray an d radio emission) to 



the light curves of the eight MSPs recently discovered by Fermi-LAT ( lAtwood et al.l 12009 
Abdo et al.ll2009dl ) in the right panels of Figures [TBI through [HH 



The Fermi-LAT light curves were produced by phase-folding LAT photons with energies 
above 100 MeV, recorded between 30 June 2008 and 15 March 2009. In order to reduce the 
contamination of the gamma-ray signal by the Galactic and extragalactic diffuse emission 
or nearby sources, and thereby maximize the signal-to-noise ratio, photons were selected in 
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narrow regions of interest, with radii of 0.5° to 1° around the pulsar locations. The gamma- 
ray light curves seen by the LAT are shown in the right panels of Figures [16] through [T9~] 
along with the radio profiles providing the absolute phase alignment. As the models predict 
different radio-to-gamma lags 5, the phase alignment is crucial. The horizontal dashed lines 
indicate the background level estimated from a ring surrounding the pulsar. 

It is important to note that the LAT angular resolution depends o n the photon energy : 
the 68% containment radius is 3.5° at 100 MeV, and 0.6° at 1 GeV (see Atwood et al.ll2009h . 
A consequence of the narrowly-chosen regions of interest is that a significant fraction of 
low-energy photons emitted by the pulsars are rejected. Therefore, the light curves shown 
in Figures [16] through [19] are biased toward energies above 1 GeV and may not reflect the 
actual profile shape obtained using larger regions of interest. 

As the Fermi mission continues, increased photon counts will allow the study of light 
curve shape as a function of energy in more detail. In fact, updated gamma-ray profiles 
based on energy-depen dent angular cut s do not differ fundamentally from what is seen in 
Figures Ql through QJ jGuillemotl l2009h . 



4. Results 



Table [D summ arizes some of the properties of the MSPs discovered by Fermi-L AT 
(lAbdo et al.l l2009dl ). All distances come from parallax measurements, except for those 



of PSR J02 



8+4232 and PSR J1614-2230 which are based on the NE2001 model (see 



effect flShklovskiilll970f ). 



Abdo et al.l (2009d N ) for references). The values of P have been corrected for the Shklovskii 



The radio beam may be quite large in the case of MSPs. Figure [6] shows examples of 
phaseplots of the radio conal beam for a = 70°. The top panel is for P = 2 ms, and the 
bottom one for P = 5 ms. The conal beam's total size and annular width become increasingly 
larger for shor ter periods, scaling as p~ - 35 . The notch, a feature of the retarded magnetic 
field solution (IDyks et al.ll2004al ). is apparent as well as increased intensity for the leading 
part, which is due to bunching of the B- field lines around the notch. (In the online version, 
plots are shown in color.) 

Differences of the geometric TPC, OG, PC, and also the PSPC models are graphically 
presented in Figures [7] and [8] (the geometric PC models do not provide particularly good fits 
to the observed light curves, and we will therefore not concentrate on their detailed properties 
in what follows). Figure [7] shows example phaseplots for TPC (top panel) and OG (bottom 
panel) models for a = 70°. For emission tangent to trailing field lines, relativistic effects of 
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aberration and time-of-flight delays cause phase shifts that nearly cancel those due to the 
curvature of the B-field, leading to accumulation of emission around narrow phase bands. 
This yields caustic structures around ~ 0.0 — 0.1 and ~ 0.4 — 0.6 in phase. (The observer 
phase is defined to be zero where the observer crosses the meridional plane which contains 
both ft and the magnetic dipole axis fi.) Emission is assumed to be symmetric for both 
magnetic poles. In the OG model, no emission originates below the null charge surface (note 
that the null charge surface is at £ = 90° in these plots), so that an observer can only see 
emission from one magnetic pole, in contrast to the SG / TPC models where an observer 
sees emission from both poles. Therefore, by comparing the two panels of Figured one can 
infer which part of the caustics originate at low emission altitudes (present only in TPC 
models), and which part at high altitudes (present in both TPC and OG models). The dark 
circular structures at phase and 0.5 in the TPC-case are the PC surfaces from opposite 
poles. They are significantly larger for MSPs than for younger pulsars, since their size scales 
with P" 1 / 2 . 

Figure [8] shows example phaseplots for the constant-emissivity PC case (top panel), and 
a PSPC model (bottom panel) including the full GR E-field, for a = 70°. The difference in 
shape of the emission regions associated with the magnetic axes in the latter case reflects 
the dependencies of the E-field on spatial parameters. The emission regions are also much 
smaller (implying correspondingly smaller gamma-ray peak widths) , as the E-field decreases 
with altitude before reaching a constant value (Figure [2]). 

In order to fit the Fermi-LAT and radio light curves (and to compare different model 
predictions; see Section [3]), we generated a large number of light curves for each of the 
different pulsar models, and for nearly the full range of inclination and observer angles 
(a = ( = 5° — 90°, in 5° intervals); also for P = 2, 3, and 5 ms, and for different gap 
widths (Section 12. 1 1) . (Although the phaseplots are usually very similar for different P in 
the case of younger pulsars, the PC size is significantly larger in the MSP case, and may 
impact light curves derived from the phaseplots.) Some example light curves generated using 
different phaseplots are shown in Figure [9] through [15] (see Table |2] for explanation of the 
model abbreviations used). Each panel shows the light curves (black: gamma-ray, gray / 
magenta: radio) corresponding to different (a, ^-combinations, with the normalized phase 
<f) = — 1 in each case. Note that all profiles have been renormalized, since we were primarily 
interested in pulse shape (and radio-to-gamma phase lag). This has the effect of boosting 
low- level emission, leading to noisy profiles in some cases (e.g., the first column of Figure [9]). 
Details as to the model, and chosen period P, are given in the captions of these Figures. 

The PSPC (and PC) model have mostly single-peaked gamma-ray profiles which are 
roughly in phase with the radio (when there is only a single radio peak), and the profiles 



-16- 



become larger when P decreases (especially the radio). The radio profile may exhibit zero, 
one or two peaks, depending on where the observer's line-of-sight intersects with the radio 
cone. In significantly off-beam geometries (large impact angle j3 = ( — a), one therefore only 
sees gamma-ray radiation (i.e., missing the radio cone), in accordance with expectations 
that gamma-ray beams are larger than their radio counterpa rts. This is the sta ndard way 



of explaining the phenomenon of 'radio-quiet' pulsars (e.g., lAbdo et al.l l2009d ). Double- 
peaked radio profiles occur for both large a and (. However, for the PSPC gamma-ray 
model, double-peaked profiles occur only for large £, because the E\\ -dependence on pc and 
r) limits emission to favorably-curved field lines at high a. Therefore, an observer mostly 
sees emission from only one pole in this case, similar to the OG model. 

Both OG and TPC models have a preponderance of double-peaked light curves at similar 
phases (see especially the lower right corners of Figures [TT1 through [T5"j) . OG models do not 
exist at all angle combinations, while TPC models do (due to emission occuring below the 
null charge surface). It is interesting to note that one may find sharp, solitary peaks for 
some regions in phase space in OG models, while the corresponding TPC-peaks usually have 
additional low-level features (e.g., compare the TPC and OG profiles at (a, C) = (30°, 60°)). 
O ur models follow the inve r se tre nd of peak separation vs. radio-to-gamma lag, first noticed 



by iRomani fc Yadigaroglul (119951 ). Our profile pulse widt h is proportional to w, because we 



assume that emission fills the full gap, unlike the case of IWatters et al.l (120091 ) 



We chose best-fit light curves from the various models to match the MSP gamma-ray 
and radio data by eye, using plots such as those in Figures 191 through [TBI However, statistical 
uncertainties in the data may complicate unique matching of predicted and observed pro- 
files. In addition, the model light curves usually do not radically change for a ~ 5°-change 
in a or (, making our obtained fits somewhat subjective. The left panels of Figures [T6l 
through [TPJ show phaseplots associated with the best light curve fits obtained for all eight 
MSPs (with horizontal lines indicating constant-^ slices). In each case, the upper subpanel 
indicates a TPC model, and the lower one an OG model, except for PSR J1744— 1134 and 
PSR J2124— 3358, where the left panels are for PSPC models. We did not find any satisfac- 
tory fits from the geometric PC models, and TPC and OG models with w — 0. In addition, 
the radio model fits the data quite well overall, except for the case of PSR J0218+4232, 
which seem to require a wider cone beam. 

In the right panels of Figures [16] to [19], we show the observed gamma-ray and radio light 
curves, along with model fits. (We normalized the data to unity. Next, we normalized the 
model light curves to unity minus the background level. We lastly added this background 
to the latter.) Three MSPs (PSR J0030+0451, PSR J0218+4232, and PSR J1614-2230) 
have double-peaked light curves, indicating the presence of screening electron-positron pairs 
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(which are necessary to form the TPC or OG emitting structure). In six cases, the gamma- 
ray light curve lags the radio. Two MSPs, PSR J0030+0451 and PSR J1614-2230, have 
a relative phase lag 5 ~ 0.2 (distinct from the function S'(rj) used earlier), and four, 
PSR J0218+4232, PSR J0437-4715, PSR J1613-0200, and PSR J0751+1807, have 5 ~ 
0.45. These MSPs are well fit by TPC and OG models. The remaining two MSPs (PSR J1744-1134 
and PSR J2124— 3358) have 5 ~ 0.85, which means that the radio lags the gamma-ray curves 
by 0.15 in phase. These two cases are exclusively fit by the PSPC model, where the gamma 
and radio emission come from the same magnetic pole, and originate well above the stellar 
surface. In the PSPC (and PC) model, the radio peak lags the gamma-ray peak, because 
the gamma-ray emission originates from all open field lines, appearing at earlier phases and 
washing out the caustic peaks. For PSR J1614— 2230, the radio profile was measured at 
1.5 GHz, and for PSR J0437-4715, at 3 GHz (although for the modeling we o nly use fre- 



quenci es 1.4 GHz and 3 GHz). All other radio profiles were observed at 1.4 GHz (lAbdo et al. 



2009dj). Our best-fit model light curves allow us to infer values for a and ( for each MSP. 
These are summarized in Table [3] (labeled with subscripts 'TPC, 'OG', and 'PSPC'), and 
compared with values obtained from radio polarimetric measurements (labeled with sub- 
scripts 'radio'). The latter inferred values are typically very difficult to obtain for MSPs due 
to the flatness of the position angle curve, and scatter of data. They are therefore generally 
quite uncertain. 

We lastly calculated fft(a,() using Eq. (134")) for each of the different models, and for 
different periods. Results are shown in Figures [2D] through [23] The 'pinpoints' of more 
intense color which are sometimes visible is an artifact of our limited resolution of 5° for a 
and (. Note that very low- level emission at large impact angles may produce excessively large 
fn factors, even for cases where the pulsar is not expected to be visible. For representational 
purposes, we set fn = when it exceeds the value of 4. We also calculated values for f n for 
our best-fit models, and summarized them in Table Although is a function of a, (, and 
P, it is typically of order unity for the best-fit geometries we consider here. The OG model 
ty pically predicts lower values than the TPC model. This is consistent with the findings 



of IWatters et all (120091 ) . Although there are small differences w hen performing a det ailed 



comparison of our results for TPC and OG models with those of IWatters et al.l (120091 ). our 
functional dependence of fn(a, £) qualitatively resembles their results, and we obtain similar 
values of fn(a, () (keeping in mind that we are modeling MSPs, while they studied younger 
pulsars). Our results for the PSPC model however differ markedly from their PC model 
results, due to the fundamental physical difference of magnetospheric structure for MSPs 
and younger pulsars (i.e., unscreened vs. screened pulsar magnetospheres) . 
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Discussion and Conclusions 



We presented results from 3D emission modeling of gamma-ray and radio radiation in 
the framework of geometric PC, OG, and TPC pulsar models, and also for the full-radiation 
PSPC model. We have applied our results to recent measurements of newly-discovered 
MSPs by Fermi-LAT . In this sense, we present results complementary to those obtained by 
Watters et al.l (120091 ) for young pulsars. 



Previously, it was believed th at most MSPs should have unscreened magnetospheres 
(IHarding. Usov. fc Muslimovl 120051 ). as they lie below the predicted CR pair death line on 
the PP-diagram. It was expected that such pair-starved MSPs should h ave single gamma- 
ray pulses roughly in phase with the radio (jVenter &: De Jagerl l2005al ) . From Figure [16] 
and [TBI we see the surprising fact that there are indeed MSPs that have double-peaked light 
curves well fit by TPC / OG models, as are many of the young gamma-ray pulsars. This is 
interpreted as indicating the operation of a magnetic pair formation mechanism, and copious 
production of pairs to set up the required emitting gap structure. 

New ways of creating pairs in low-P rot pulsars will have to be found to explain this 
phenomenon. PSR J0030+0451 illustrates this point very well in that it has the lowest E TOt 
of the MSP sam ple (3.5 x 10 33 erg s~ l ), therefore lyin g significantly below the calculated CR 



death line (e.g., IHarding. Muslimov. fc Zhang] 120021 ). and yet exhibits the sharpest double 



peaks of the current population, implying emission originating in very thin TPC / OG 
gaps. The problem may be alleviated somewhat by increasing the stellar compactness k' 
(larger mass or small er radius), motivated by recent measu rements of large MSP masses 



(up to ~ 1.7M Q ; see lYerbiest et al.l 120081 ; iFreire et al.ll2009l . and references therein). This 



will boost the GR E-fields, and enhance pair creation probability. Another way to do this 
would be to increase the magnetic field. B-fields that are larger than those usually inferred 
using the dipole spin-down model (and having smaller curvature radii) may be present when 
there are multipolar B-components near the surface (or an offset-dipole geometry). In fact, 
offset dipoles have been suggested i n modeling the X-ray light curves of MSPs J0437— 4715 
and J0030+0451 faogdanov et all l2007i : bogdanov fc Grindlavl l2009h . However, detailed 
investigation of such a scenario and its implications for pair cascades is necessary to place 
this speculation on sure footing. Another possible origin for higher surface fields is the 
movement of mag netic poles toward the spin axis during the spin-up phase of an MSP 
(ILamb et al.l 120081 ) . It has been argued that during the spin-up to millisecond periods, the 
inward motion of the neutron star superfluid vorti ces produces a st rain on the crust, causing 
the magnetic poles to drift toward the spin-axis (jRudermanlll99ll ). If the two poles are in 
the same hemisphere prior to spin-up, the poles drift toward each other, producing a nearly 
orthogonal rotator having the same dipole moment but a surface field that can be orders of 
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magnitude higher (IChen fc Rudermanl 119931 ). 



We find that there is exclusive differentiation between the TPC / OG models on the one 
hand, and the PSPC model on the other hand. Six MSPs have gamma-ray light curves which 
lag the radio and are explained using TPC or OG fits, but not PSPC fits. For the remaining 
two MSPs, the radio light curves slightly lag the gamma-ray light curves, and these are fit by 
the PSPC model (and not by the TPC / OG models). It therefore seems that there are two 
subclasses emerging within the current gamma-ray MSP sample, and it is not obvious which 
pulsar characteristics provide a means to predict subclass membership. From our model 
light curve fitting, we furthermore find (a, Q values which are in reasonable agreement with 
values inferred from MSP polarization measurements. Although we find good PSPC fits 
for the last two MSPs, we caution that the E-field is only approximately known (e.g., it 
follows from a local electrodynamical model based on a GR dipolar B-field). Future models 
which take global current flow patterns into account, along with more sophisticated B-field 
structure, may produce more realistic solutions for the E-field. 

Our ability to discriminate between different classes of models derives from the fact that 
we produced both the gamma-ray and radio curves within the same model. We could then 
use the shape and relative radio-to-gamma phase lag provided by the data to obtain the 
best-fit model type for each MSP. The data also enabled us to conclude that the emission, 
in all models considered, must come from the outer magnetosphere. This has now be en 
observed to be true for the bulk of the gamma-ray pulsar population (lAbdo et al.ll2009el ). 



In the case of PSR J0437-4715 and PSR J06 13-0200, we find that the TPC model 
predicts a significant precursor to the main gamma-ray peak, while the OG model predicts no 
such low-level emission. With more statistics, this effect may possibly become a discriminator 
between the TPC and OG models. (We assumed that the TPC emission region starts at 
r em = R when creating our plots. However, the relative intensity of the precursor and low- 
level emission predicted by the TPC model may be reduced by limiting the emission region's 
extension, i.e. only collecting photons above a certain radius r em > -R m i n > R-) 

We calculated the flux correction factor in the context of the different models, and found 
that /n ~ 1. These values imply a wide beaming angle , and derives from the fact that we 
obtain best fits for large impact angles. IVenterl (120081 ) previously found A avg ~ 10 — 30 
(i.e. f a ~ 0.8 - 2.4), and A max ~ 300 (/^ ax ~ 24) for the PSPF model. Now, we find fa 
~ 0.5 — 2, and fa*** ~ 4 . These result s are roughly consistent, with the differences stemming 
from the following: (i) IVenterl (120081 ) used energy flux ratios to calculate A, while we are 



using photon flux ratios to calculate fa, ass uming that th e photon and energy fluxes have 
similar distributions across (£, </>)-space; (ii) IVenterl (120081 ) only used E, 



and E^ for the 



E-field, while we now also include the high-altitude solution (E^) for the PSPF case. This 
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leads to more intense high-altitude emission, and therefore smaller values of /q for off-beam 
emission. 

We lastly remark that the larger radio beam widths of MSPs compared to those of 
canonical pulsars should lead one to expect relatively few radio-quiet MSPs. 

The spectacular data from Fermi-LAT hold the promise of phase-resolved spectroscopy, 
at least for the brightest pulsars, and will challenge existing pulsar models to reproduce such 
unprecedented detail. Future work therefore includes using full acceleration and radiation 
models to study gamma-ray spectra, luminosities, and light curves, in order to constrain 
fundamental electrodynamical quantities, and possibly providing the opportunity of probing 
the emission geometry and B-field structure more deeply. Improved understandi ng of pulsar 



mode ls will also feed back into more accurate population synthesis models (e.g.. IStory et al. 



20071 ). In addition, we hope to obtain better understanding of important quantities such as 



MSP e fficiencies, and whe ther this quantity is similar for Galactic-Field and globular-cluster 



MSPs flAbdo et al.ll2009d l. 
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Table 1. Parameters of MSPs discovered by Fermi-LAT ( lAbdo et al 



2009d) 



Name 


P 


P 


Distance 


Age 


E ro t 


B 




(ms) 


(io- 20 ) 


(kpc) 


(10 9 yr) 


(10 33 erg s -1 ) 


(10 8 G) 


J0030+0451 


4.865 


1.01 


0.300 ± 0.090 


7.63 


3.47 


2.04 


J0218+4232 


2.323 


7.79 


2.70 ± 0.60 


0.47 


245 


4.31 


J0437-4715 


5.757 


1.39 


0.156 ± 0.002 


6.55 


2.88 


2.87 


J0613-0200 


3.061 


0.915 


0.48 ± 0.14 


5.31 


12.6 


1.69 


J0751+1807 


3.479 


0.755 


0.62 ± 0.31 


7.30 


7.08 


1.64 


J1614-2230 


3.151 


0.397 


1.30 ± 0.25 


12.6 


5.01 


1.13 


J1744-1134 


4.075 


0.682 


0.470 ± 0.090 


9.47 


3.98 


1.69 


J2124-3358 


4.931 


1.21 


0.25 ± 0.13 


6.47 


3.98 


2.47 
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Table 2. MSP Model Descriptions 



Abbreviation 


Tovc 


w 




Azimuthal bins 


Description 


TPC1 


[0.90,1.00] 


0.10 


0.005 


180 


Geometric TPC Model 


TPC2 


[0.95,1.00] 


0.05 


0.005 


180 


Geometric TPC Model 


TPC3 


[0.80, 1.00] 


0.20 


0.005 


180 


Geometric TPC Model 


TPC4 


[0.60, 1.00] 


0.40 


0.005 


180 


Geometric TPC Model 


TPC5 


[1.00,1.00] 


0.00 


0.005 


180 


Geometric TPC Model 


OG1 


[0.95,1.00] 


0.05 


0.005 


180 


Geometric OG Model 


OG2 


[0.90,0.90] 


0.00 


0.005 


180 


Geometric OG Model 


OG3 


[1.00,1.00] 


0.00 


0.005 


180 


Geometric OG Model 


PCI 


[0.00, 1.00] 


1.00 


0.005 


180 


Geometric PC Model 


PC2 


[0.00, 1.00] 


1.00 


0.005 


180 


Radiation PSPC Model 
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Table 3. Model fits for a, (, and fn(at, (, P) 



Name 


QTPC 


Ctpc 


QOG 


Cog 


«PSPC 


Cpspc 


^radio 


Cradio 


Ref. 


/n,TPC 


/n,OG 


/n.pspc 




(°) 


(°) 


(°) 


(°) 


(°) 


(°) 


(°) 


(°) 










J0030+0451 


70 


80 


80 


70 






~ 62 


~ 72 


1 


1.04 


0.90 




J0218+4232 


(50 


60 


50 


70 






~ 8 


~ 90 


2 


1.06 


0.63 




J0437-4715 


30 


60 


30 


60 






20 - 35 


16 - 20 


3,4 


1.23 


1.82 




J0613-0200 


30 


60 


30 


60 






small j3 




5 


1.19 


1.76 




J0751+1807 


50 


50 


50 


50 












0.80 


0.65 




J1614-2230 


10 


80 


10 


80 












0.83 


0.64 




J1744-1134 










50 


80 












1.19 


J2124-3358 










40 


80 


20 - 60 (48) 


27 - 80 (67) 


6 






1.29 



References. — (1) iLommen et alj j2OO0h; ( 21 IStairs et al.l l|l999h : (3) iManchester fc Johnston! l|l995h : (4) iGil fc Krawczvkl jl997h . (5) 
IXilouris et al] lll998l l; C61 IManchester fc Hani j2004l 
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Fig. 1. — The PP-diagram, indicating contours of constant E m t (dashed lines) and rotational 
age (solid lines), as well as pulsars from the ATNF Catalog ([Manchester et al.ll2005l ). We used 
values of P > corrected for the Shklovskii effect (jShklovskiil 1 19 70l ) . an d removed pulsars i n 
globular clusters. The squares are the 8 newly-discovered Fermi MSPs (lAbdo et al.ll2009dl ). 
All except PSR J021 8 +4232 lie below the ICS d eathline, and all eight lie below the CR 
deathline (modeled by lHarding fc Muslimovl |2002| ) . 
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Examples of the final E-field we obtain after matching e9~' through for 
irameters: —Eu vs. log 10 of the height above the PC, normalized by the PC 
= (QH 3 /c)^ 2 . These plots were obtained for P = 5.75 x 10~ 3 s, P = lO" 20 , 
0. in panel ya.), we cnuse u. = z,u , c; = u.o, yj pc = 45°, in 

200°. 



R = 10 6 cm, and M 
panel (b), a = 35°, £ 



1AM, 
0.7, 



'pe- 



rn panel (a), we chose a = 20°, £ = 0.3, <f) pc 
= 150°, and in panel (c), a = 80°, £ = 0.8, <f> 



pc 



In the last panel, the final —E\\ is negative, so no solution of rj c is obtained. In each case, 
we label through (thick solid lines), indicate potential solutions (which vary with 
r) c ) by thin gray (cyan) lines, and the final solution by thick (red) dashed lines. Also, we 



indicate r/b where we match and E^\ 



(2) (3) 

and 77 c where we match E^ and E^ , by thin 

,(3) 



vertical dashed lines. (Although E^' does slightly vary with r) c: we only indicate the dis- 
solution corresponding to the r) c found for the final solution. For panel (c), we show a typical 
E, i -solution.) 




Fig. 3. — Contour plots of our solutions of 7] c for P = 5 ms, and for a = 10°, 20°, ...,90°; 
£ is the radial and pc the azimuthal coordinate in each case. The magnetic dipole axis [i 
is situated at the origin, pointing outward normal to the plane of the page, in each case. 
The rotation axis f2 is in the direction of pc = 0, while the leading (trailing) edge of the 
pulse profile originates on B-field lines with footpoints around (fi pc ~ 90° (</> pc ~ 270°). The 
^-solutions get progressively smaller for <ft pc ~ 180°, and for large a, until no solution is 
found which satisfies our solution matching criteria (denoted by zero values or no values 
at all on the plots above). We ignore the emission from those particular field lines. We 
expect the ^-distribution to reflect the symmetry of the cos0 pc function which is found in 
E\\; the small irregularities stem from the fact that we used interpolation on a non-uniform 
(£> 0pc)-grid when preparing the contour plots. (See online version for color plots.) 
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Fig. 4. — The log 10 of gain (acceleration) rate 7 gain (solid line), loss rate 7i oss (dash-dotted 
line), curvature radius p c (dash-dot-dotted line), and the Lorentz factor 7 (short-dashed line) 
as function of normalized radial distance r\. We used pc = 360°, £ = 0.7, a = 40°, and 
P = 5 ms in this plot. 
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Fig. 5. — The log 10 of cumulative energy gain (J =1 rf7 g ainj solid line), cumulative energy- 
losses (JZLi ^7ios S ; dash-dotted line), Lorentz factor 7 (short-dashed line), and the sum of the 
cumulative losses and acquired particle energy (P =1 d^\ oss + 7(77) — 7 ; thin dashed gray / 
cyan line) in units of m e c 2 vs. rj. The latter sum and the cumulative gain coincide (within 
~ 0.3% for the //-range shown), pointing to energy balance, i.e. electric potential energy 
being converted into gamma-radiation and particle kinetic energy. We used pc = 360°, 
£ = 0.7, a = 40°, and P = 5 ms for this plot. 
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Fig. 6. — Example phaseplots of the radio conal beam, for a = 70° and at a frequency of 
1.4 GHz. Panel (a) is for P = 2 ms, and panel (b) for P = 5 ms. Beam and annulus widths 
become increasingly larger for shorter periods. Increased intensity for the leading part is 
due to bunching of the B-field lines around the notch. (Note that in this and following 
phaseplots, the color scales are not the same for the different panels, but are chosen to show 
the most details for each case.) 
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Fig. 7. — Example phaseplots for the TPC2 and OG1 models (panel (a) and (b) respectively), 
for a = 70° and P = 5 ms. In contrast to the TPC models, no emission originates below the 
null charge surface in the OG model, in which case an observer can only see emission from 
one magnetic pole. 
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Fig. 9. — Sample light curves (black: gamma-ray; gray / magenta: radio at 1.4 GHz) for 
the PC2 model, with P = 2 ms. The observer angle ( changes along the columns, and the 
inclination angle a along the rows. All pulse shape maxima are normalized to unity, and the 
phase range goes from = — 1 in each case (and similar for subsequent figures). 
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Fig. 11. — Sample light curves for a TPC1 model with P = 2 ms. 
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Fig. 12. — Sample light curves for a TPC2 model with P = 3 ms. 
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Fig. 13. — Sample light curves for a TPC2 model with P = 5 ms. 
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Fig. 14. — Sample light curves for an 0G1 model with P = 2 ms. 
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Fig. 15. — Sample light curves for an OG1 model with P = 5 ms. 
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Fig. 16. — Gamma-ray phaseplots (left panels) and observed and fitted gamma-ray and 
radio light curves (right panels) for PSR J0030+0451 (panel [a]-[d], P = 5 ms) and 
PSR J0218+4232 (panel [e]-[h], P = 2 ms). Panel (a) is for a TPC1 model with 
(a, C) = (70°, 80°), (b) for an OG1 model with (a, () = (80°, 70°), (e) for a TPC1 model with 
(a,C) = (60°, 60°), and (f) for an OG1 model with (a, C) = (50°, 70°). For the gamma-ray 
light c urves (panel [c] and [g]), the histograms represent the Fermi-LAT data (lAbdo et al. 
2009dJ), the horizontal dashed line the estimated background level, the dashed (online: ma- 
genta) lines are TPC fits, and dash-dotted (online: green) lines are OG fits (see Table [3]) 
For the radio light curves (panel [d] and [h]) the solid (blue) line represents the radio data 
while the dashed (magenta) and dash-dotted (green) lines correspond to the same (a, ()■ 
combinations as those of the respective TPC and OG fits. 
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Fig. 17.— Same as Figure Ml but for PSR J0437-4715 (panel [a]-[d], P = 5 ms) and 
PSR J0613-0200 (panel [e]-[h], P = 3 ms). Panel (a) is for a TPC2 model with (a, () = 
(30°, 60°), (b) for an OG1 model with (a, C) = (30°, 60°), (e) for a TPC2 model with (a, () = 
(30°, 60°), and (f) for an OG1 model with (a,() — (30°, 60°). For the cases where we use 
the same (a, (^-combination for both the TPC and OG fits, we only have a single radio light 
curve fit. 
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Fig. 18.— Same as Figure M, but for PSR J0751+1807 (panel [a]-[d], P = 3 ms) and 
PSR J1614-2230 (panel [e]-[h], P = 3 ms). Panel (a) is for a TPC2 model with (a, () = 
(50°, 50°), (b) for an OG1 model with (a, () = (50°, 50°), (e) for a TPC2 model with (a, () = 
(40°, 80°), and (f) for an OG1 model with (a, () = (40°, 80°). 
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Fig. 19.— Similar to Figure [T6l but for PSR J1744-1134 (panel [a]-[c], P = 5 ms) and 
PSR J2124-3358 (panel [d]-[f], P = 5 ms). Panel (a) is for PC2 model with (a, () = (50°, 80°), 
and (d) for a PC2 model with (a, () = (40°, 80°). In panels (b) and (e), the dashed (magenta) 
lines signify PC2 model fits. 
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Fig. 20. — The flux correction factor fa(a, £, P) vs. a and ( for a TPC2 model, with panel (a) 
and (b) for P = 2 ms and P = 5 ms, respectively. The same color scale is used throughout, 
and values of /q > 4 were set to zero. 
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Fig. 21. — The flux correction factor fa(a, (, P) vs. a and ( for an 0G1 model, with panel (a) 
and (b) for P = 2 ms and P = 5 ms, respectively. 
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Fig. 22. — The flux correction factor fn(a, (, P) vs. a and £ for a PCI model, with panel (a) 
and (b) for P = 2 ms and P = 5 ms, respectively. 




Fig. 23. — The flux correction factor fn(a, (, P) vs. a and ( for a PC2 model, with panel (a) 
and (b) for P = 2 ms and P = 5 ms, respectively. 



